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Abstract 

A fully conservative numerical method for the computation of steady 
invlscid supersonic flow about general conical bodies at incidence is 
described. The procedure utilizes the potential approximation and Implements 
a body conforming mesh generator. The conical potential is assumed to have 
its best linear variation Inside each mesh cell; a secondary interlocking cell 
system is used to establish the flux balance required to conserve mass. In 
the supersonic regions the scheme is desymmetrizled by adding artificial 
viscosity in conservation form. The algorithm is nearly an order of a 
magnitude faster than present Euler methods and predicts known results 
accurately and qualitative features such as nodal point lift off correctly. 
Results are compared with those of other investigators. 
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Introduction 


Conical flows are one of the simpliest type of inviscid flows that 
have the basic features of a three-dimensional flow. A flow field is 
classified as conical when all the physical properties, viz., the 
pressure, density, velocity, and entropy, remain constant along every 
straight line through a given point called the apex. Conical flows 
occur for example, around finite cones in supersonic flows because of 
the law of forbidden signals. The topological features of conical flows 
can be easily understood by studying the cross-flow streamlines; that is 
the traces of the conical stream surface's intersection with a sphere as 
sketched in Figure 1. The cross-flow streamlines will have critical 
points where the cross-flow velocities vanish. For a special class of 
critical points one can derive rules for the number of these points 
using Poincare indices. Thus, for example, irrotational conical flows 
must have an equal number of saddle points and nodes. At a nodal point 
the entropy, density and radial velocity are multivalued. At high 
angles of attack conical streamline patterns exhibit certain global 
changes such as the lift-off of the leeward node and, perhaps, the 
appearance of spiral nodes. In addition, the cross flow may become 
supersonic as it expands about the leading edge, leading to an embedded 
supersonic cross-flow region terminated by a shock wave (see Figure 1). 
It is interesting to note that the perturbation of shock free flows to 
see if neighboring solutions exist in the classical sense of Morawetz^ 
remains to be studied for conical flows. 

The isentropic assumption retains all of the topological features 
of these flows except spiral nodes and should provide an adequate 
approximation to the quantitative flow features if the Mach number 
normal to any shock is less than about 1.4. This approximation greatly 



simplifies the computations because the governing equation is scalar and 
the possibility of multiple values at a nodal point (the vortical 
singularity) is eliminated. 

First major success in computing nonlinear irrotational conical flows was 

O 

reported by Grossman'^, who devised a quasilinear finite difference method to 
this problem. An alternative approach is to extend Jameson's^ finite 
difference algorithm for transonic full potential equation in the euclidean 
three space as a marching scheme to treat supersonic potential flows. This 
extension has been developed by Shankar^ using Steger's^ density linearisation 
technique . 

In this paper the theory of irrotational conical flows is described in a 
general coordinate system defined on a unit sphere. A finite area method is 
described that represents an extension of the conventional finite volume 
methods to vector fields defined on a curved surface. It is then used to 

on in 

compute various conical flow fields studied by other investigators * 

Simlarlty between the hlghest-order terms of the partial differential equation 
describing conical flows and that describing plane transonic flows has been 
exploited to devise a suitable artifical viscosity to implement the 
(mathematical) entropy condition^ ^ as well as to construct a stable Iteration 


scheme . 
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Cross-flow velocity field 

It IS essential in the application of the finite area method to 
frame the equations in an invariant coordinate system. We have done 
this by first projecting the Euler equations for a general 
three-dimensional flow onto a sphere of radius r, and then scaling them 
to obtain the description on the unit sphere. We first note that the 
main stream velocity components, Q'' , and their projection on the 
sphere of radius r, are related by 


V“ = B“Q^’ 


where 


B 


i 

a 



as 


a 


1,2 and i 


1,2,3 


A 4C. . 

are the projection factors (or the tangents). Here the X'' are 
coordinates in Euclidian 3-space and the 5 “ are parametric coordinates 
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on the surface of a sphere. We could also decompose the mainstream 
velocity as 


qT = B^V“ + QdN’ 

Ct K 

where Qp is the radial velocity and the normal to the spherical 
surface. 

Now, with the conical assumption in mind, we introduce a radial 
scaling to reduce the variables to those on a unit sphere, viz., 

V 

V“ = rV“ , 

such that the magnitude of the scaled cross-flow velocity V“V = V“V = 

a a c 

is independent of r. The total velocity is thus 

Once the coordinates 3“, the metric g^^p, and the velocity vector 
Y“, are defined on a portion of the surface of the unit sphere, we may 
use elementary results from Riemann geometry to develop the governing 
equations. For potential flow we need to use only the continuity and 
the energy equations. 

The continutiy equation for the general three-dimensional inviscid 
flow is 


6p/G q" 

sx^ 


0 
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where /G is determinant of the metric tensor of the space , and p 
the gas density. For conical flows this implies that on a unit sphere 

+ 2p/g q„ = 0 

as” “ 

Here /g is the determinant of the metric tensor of the parametric space 

ff 13 

s on the surface of unit sphere. If the irrotational assumption 
is made, the velocity Q'' will have a potential such that for 

conical flows 


(|)(XM = rF(3“) 

where F(s“) may be called the conical velocity potential since 


V = — , Qn = F(s“) 

a K 

0^ 


We also have the energy equation 

pY-l = 1 + m2 (1-Q2) , 


( 1 ) 


where 


q 2 = v“v + f2 

a 


and M is the freestream Mach number. Thus, we must solve the equation 


5p/g r 


+ 2p/g F = 0 


( 2 ) 




V“ = g“P ^ 

asP 


where 


(3) 
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and the density is given by the energy equation (1). 

Substituting (1) and (3) into (2) and performing the differentiation, we 
find the quasi! inear form of the governing partial differentia! 
equation: 

p/g [(9“** - v^,p + ( 2 - F] = 0 Ml 

1 9 

where "II" means surface covariant differentiation^'^ and, a, the sound 
speed and is given by, 

• pl-l (5) 

03 ^ 

Equation (4) changes its type when 


12 

(g 



11 


22 


(g (g ) 


= 0 


or, noting g > 0 always, when 


1 

g 



- 1 ] = 0 


Here we use the notation 


V“ = (ll), V = (“) and s“ = th 

V a V p 

Thus, the equation is hyperbolic, parabolic, or elliptic depending on 
whether the cross-flow Mach number M,- is the greater, equal, or less 
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than one. We could also derive this result by choosing a local 
coordinate system aligned with the cross-flow streamlines. 

If we define streamwise and normal coordinates s,n such that 


a an 


= 0 and e „V 


a 53^ 


ap as 


= 0 . 


with 


_ ^ - a 

^ap as as ^ 


Oi Oi 


ap 5n an 


= 1 


1 2 

Here is the surface permutation s 3 rmbol 


Then we find 


as“ _ v“ a^ _ 1 


as 


qc ’ an 


^9 q 


^ and ^ = 


an _ 1 u 


an 


/g q. 


Note also that in these coordinates 


aF _ aF 


bi 


as .„a as 
a^ 


= q. 


V r 

a 


and 


a_F ^ aF ( _ 1 


an 


a^ 


_ 1_ ) + ^ (_i 1_) 
/g q^, an >^9 q^ 


= 0 


Similarily, we also obtain the relations 


a^F _ V"V 




a^F 


as" 


q2 as^as' 


+ ... 


and 


a^F ^ (g«p _ ^ a^F 

an2 


*^c as“a3^ 



Using these relations we may write the governing partial differential 
equation as 


- [(q2-a2).^^ -a2^]+... =0 

2 c 


6S' 


dn‘ 


We also note the following structure of these equation: 
the expression for 5^F/ds^ we obtain 


- 2 :^ CnV«VP 
=2 


0.*i 


2 

511 "^ 


= 0 


where 



Now if we define two functions P, Q such that 

P = (ip (U2 ^ + UV ^LL- ) 
a s^an 

and 


Q = pp (UV + V2 

a a^dri dn 

then the partial differential equation becomes 


-(P + Q) + p/g ^ + 

dn^ 


0 


Replacing 


This form is useful for explaining the introduction of a conservative 
artificial viscosity. 
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Finally, we note that the governing equation is of quasi! inear type 
and hence, that it admits shock jumps. The jump condition that 
conserves mass follows immediately from the Eq. (2): 

(^) = Mi 

^ shock iipUii 

This is obviously the jump condition that would be desired from the 
mainstream continuity equation with the conical assumption. 

Finite area method 

O 

Jameson's finite volume method for the potential equation may be 
extended to a vector field defined on a non-Eucl idi an space as long as 
we have a similar partial differential equation. In this section we 
develop the finite area method on the unit sphere. It should be 
emphasized that the derivation would be the same for a vector field on a 
general curved surface. We assume that on this curved surface a smooth 
grid is provided, as sketched in Figure 2, with the surface coordinates 
(latitude e, longitude cp) provided for each nodal point: 

e“ = (®) 

We call these primary cells. In order to implement the finite area 
method the primary cells are mapped to a unit square using a local 
bilinear transformation in the parametric space such that 

4 . . 4 . . 

= E s’e’ and cp = E 
i=l i=l 


0 
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Here i denotes the nodal values and 


S^' = 4(i + + n^'n) 


The geometrical quantities are calculated in the following manner. 
The first fundamental form in the spherical coordinate system e« is 
ds^ = sin^4^d9^ + dcj;^ , or 


sin^ci; 0 

q = ^ 

0 1 


( 6 ) 


and 


/g = sin (1; 


( 7 ) 


In a mapped coordinate system with ds^ 


g ds'^dH^, then 


and 


■^ap 


- 80 ^ 89 ^^ 

8S“ 8hP 


/g = /g J 


( 8 ) 


Here J is the Jacobian of the parametric transformation 0“ (z^), that is 


/g = sin 4. (0.(1; - 0 ();_) ( 9 ) 

We always calculate the geometric quantities at the center of the cells 
and therefore the bilinear transformation and its best linear substitute 
have the same role^\ Thus we take 
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and thus 


ci - ci - 1 X r''r 4. ^ 

s + CS + Tiri 


e_ = E 0^5’ , etc. 
^ i=l 


and at the center of the cell 


1 ^ 

0 = 7T E 0^ etc. (10) 

^ i = l 

Equations (6) - (10) define geometric quantities at the center of 
primary cells. The flov/ quantities are also defined at the center of 
the cells. The potential is of course defined only at the nodal 
points. The flow quantities may be calculated as follows: Let f be the 

disturbance to the freestream potential f„ due to the body, i.e., F = 
f„ + f. Then we assume the disturbed potential also to have the 
bilinear form 


4 

f = E 
i=l 

Because lumping is not used in the present formulation (it was not found 
to be necessary), v;e may replace by . 


4 

f = E 
i=l 


^b 


Thus, 
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4 . . 

and f = E , etc. 

^ 1=1 

The total velocity is computed from 

v« = g»p [-Z + _5f] 

5H^ 5S^ 

and 

/V j\ 

= cos a COS 4 ; + sin 0 sin 4 , sin a, 
where a is the angle of attack. 

We are now ready to implement the finite are method. We first 

introduce a secondary interlocking cell structure as shown in Figure 3 

in order to integrate the mass continuity equations. We now integrate 
the weak conservation law over the domain q so that 

/ 7a /gdC^n + / 2pF/gd^dn - 0 

Q ^ as“ Q 

Applying the surface divergence theorm to the first term, we find 

/ p/g V“n dS + /2pF/gd^dn = 0 
c “ Q 

Here is tangent to the surface and normal to the curve c • This relation 

is valid for any arbitrary ^ and therefore also valid locally for a flux 

cell. Since the flux cell faces are parallel to coordinate lines in the 
mapped plane, and using one point evaluation for each integral, we obtain 


6[pU/g] + 6[pV/g] + (2pF/g)p = 0 
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for each cell, where 6[...] denotes the net flux change across the cell 
and (..-)o average over the cell. We also need to define the flux 
quantities at points A,B,C and D of the flux cell faces. In this 
formulation we simply use a box scheme to evaluate these terms. Thus we 
obtain the approximation to Equation (2) as 

p 6_(p/gU) + p 6 (p/gV) + p p (2p/gF) = 0. 
n 5 5 h ^ T) 

where p and 6 are respectively the averaging and central difference 
operations. 

Boundary conditions 

We consider the computational domain in Figure 4. The outer 
boundary Cq is taken well outside the bow shock wave. Boundaries and 
C 2 are symmetry planes and Cb is the cone body where the normal 
velocity vanishes. 

Outerboundary 

At the outer boundary all the disturbance vanish, i.e., f, f^, are 

all zero. This is implemented in the following way, if N 2 grids are the 
rings , then 

fCl.N^) = fd.N^ + 2) and f(I,N 2 + 1) = 0 

Symnetry plane 

At the symmetry plane v/e introduce an additional grid line and 
explicitly set the reduced potential to be the same on both grid lines. 
Thus, if Nj grids are in the circumferential direction, then 


f(l,J) = f(3,J) and f(Ni - 1,J) = f(Ni + 1,J) 
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Cone surface 

On the body surface the normal velocity should be zero. If ^(s“) 

Is the cone surface, then U^9^/3X^ = 0 implies v“3^/35“= 0 since the 

body is a cone. If the body coincides with a coordinate surface, ^ for 
example, then = n = 0 and the boundary condition implies V = 0, 

i.e., the contravariant cross-flow component that does not lie on the 
body must vanish. This is implemented by considering a half flux cell 
about the cone and using flux reflection. 

Artificial viscosity 

In order to stabilize the scheme in the supersonic regions we de- 
symmetrize the scheme by upwind differencing the contribution for the 
Fss term. Also, since the higher-order terms of the partial differen- 
tial equation for conical flows are similar to that of plane transonic 
flows, if we do the upwind differencing with first-order accuracy (at 
least near the shocks), then the resulting truncation errors will look 
like the viscous term for plane flows and therefore may be expected to 
capture any shock waves and insure the entropy condition. The viscosity 
should be Introduced in the conservation form and this can be accomplished in 
the following manner. Let us consider the case when > 0. We noticed 

earlier that the terms contributing the Fgg term have a structure containing 

-(P + Q) and therefore will effectively evaluated as -P. . “ Qj . 

^ > J ^ > J 

finite area scheme. To upwind we need to replace 


in the 
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thein by -P-_, 4 - Q,- 4 i and this means we need to add 

Tfj- ' 'N.j - 

R.j. We use the switching function 


a viscosity term 
to the residual 



in the supersonic zones 
in the subsonic zones 


The method has to be appropriately modified for the other directions of 
contravariant velocity. 

Iteration scheme 

The nonlinear algebraic equations resulting from the finite area 
method are solved using a "constructed" line relaxation scheme. This 
means that v/e assume that the problem is being solved by Jameson's 
rotated difference scheme in the quasi linear form along with Jameson's 
special relaxation method^. Here we assume that the iteration process 
is equivalent to a problem of evolution in an artificial time and choose 
the explicit time dependent terms such that the problem is well posed. 

We do have certain restrictions on the direction of the sweep. We 
should not sweep against the flow inside the supersonic zone. At high 
angle of attack this condition is difficult to maintain with a ring 
relaxation inward from the bow shock wave. Thus the suitable line 
relaxations are either a circumferential sweep or a combination of the two 
(see Figure 4). One should note that the restrictions on the sweep direction 
can be easily removed by devising an approximate factorization scheme^^. We 
use line relaxation mainly to test the finite area formulation. It was found 
that a circumferential relaxation from windward to leeward is the best in most 
of the cases and is this scheme described here. Assume tht > 0 and 


consider the line relaxation scheme. 
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scheme 


AiC.j + AsiC.j - + A3(C.j - + A^lC.j - 

" "’’ij ■ ^^i-l,j+l^ 

If we now consider that R-. + T.^ are equivalent to their quasi linear 
finite difference equivalent multiplied by p/g, then, construction of 
the Jameson iterative scheme will give the following values for 

Ai » • • • jA5 : 

, . ii2 o 

subsonic 


A, = 


p/g{g^^ - \)(- -1) 

a w 


0 


supersonic 


where w is the over-relaxation factor, 


A2 = p/g(g^^ 





A3 = p/g(g^^ 


v2 

a^ 




p/g(g 


22 




and 


A^ = 


- p/g(g 


12 


uv 


yv] 


+ 


where p is the switching function. We note here that in subsonic flow, 
provided 0 ) < 2, all the coefficients are positive. Thus the 
scheme is linearly stable. We obtain further insight by looking at its 
equivalent time dependent form for each flow type; 
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Subsonic zone 




‘1 - ''c’Fss 


+ F + 
nn 


Supersonic zone 


Atr - F - (l-M^)F 1 = (1 - M^)F + F + 


In the supersonic zone we apply the condition, 


[^<«c - ^ “ 


to ensure that s is time-like in the unsteady problem as well. This 
means that 


g(M2 - 1) > 1 

To ensure that the above condition is always satisfied, especially 
near the sonic line, we further augment the first term by adding 
e(U + V) F^^ /q^ where e is as small as possible and yet sufficient 
to ensure stability. The term F^^ has to be represented by an upwind 
difference, we write this as 


/U + V. pU 





C_. 


i J-1 


)] 


This scheme must be appropriately modified when V changes sign. 
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Grid generation on a curved surface 

Suppose we are interested in generating grids for a simply- 
connected region on a curved surface. We could use the following simple 
method: First we note that the first fundamental form is in the e“ 

coordinate system is 


ds^ = g „d0“de^ 

ap 

We first transform to a new coordinate system h“ such that 

ds^ = + dn^) 

where X = X(C,n) 

This coordinate system is called the isothermal coordinate system and 
this transform maps the surface portion conformally to a plane. Then we 
define a complex variable z such that 

z = I + ir) 

and apply further conformal transformations to obtain a simple domain 
where we may generate the desired grids. Alternately, one could derive 
a numerical grid generation method for the isothermal coordinates. In 
the present problem we used the grid generation procedure that is 
commonly used in supersonic computations, that is, we obtain the 
isothermal coordinates for a unit sphere using stereographic projection 
and then use a Joukowski transformation followed by a simple shearing to 
obtain a suitable grid network. 
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Results and discussion 

Computations were made to demonstrate that the method predicts 
qualitative features of simple flows correctly and their quantitative 
aspects accurately. All the calculations were performed on three mesh 
levels start! nq with a 16 x 16 grid system. On this initial grid 150 
iterations were performed and this was followed by 100 iterations at the 
32 X 32 level and at the final level. Convergence for the last two 
grids is reliable after 25 to 50 iterations, depending principally on 
whether or not there is a body shock wave. Calculations were performed 
with uniform grids without any clustering; a typical grid is shown in 
Figure 5. The results for a circular cone of 10° half angle at 10° 
angle of attack in a freestream of Mach 2, are shown in Figures 6a and 
6b. While the results shown are for a 64 x 64 grid, excellent agreement 
for the pressure coefficient on the body and the bow shock positions was 
obtained using the 16 x 16 grid and this required fewer than 70 itera- 
tions. This coarse grid only requires few seconds of COC 7600 CPU 
time. The pressure distribution in the field is shown in Figure 7 for 
three circumferential angles. Excellent agreement with the Euler compu- 
tations of D.J. Jones ^ is again demonstrated. An example of lift off is 
given in Figure 8, where the streamline patterns for the 10° angle of 
attack case are compared with those for an angle of attack of 20°. 

Results for tv/o thin elliptic cones are shown in Figures 9 and 10. 
One ellipse has a major to minor axis ratio of approximately 6:1 and the 
other has a ratio 13:1. In the first case a comparison with the Euler 
equation calculations of Siclari^°is made. The agreement is generally 
excellent except for the extra leading edge suction which may be due in 



part to the potential approximation, and except for the post shock 
pressure. We note here that the Euler result does not show the expected 
shock foot singularity* captured by our potential calculations. The 
second example compares the Euler equation results of Siclari, the non- 
conservative potential finite difference results of Grossman^ and our 
results. Here all three methods capture to some extent, the shock foot 
singularity. The finite area method agrees well with the Euler 
results. The difference in the shock position between the conservative 
and nonconservative method is to be noted. The total computation time for 
case with body shock wave is about 40 seconds of CDC 7600 CPU time. 


'k 

The pressure gradient immediately behind the shock must be 
logarithmically infinite. 
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PROJECTION OF STREAMLINES 
ON SPHERE R = CONST- 


PARTICLE TRAJECTORY 


CON I CAL 
.STREAM SURFACE 



SHOCK WAVE 


BOW SHOCK WAVE 



CONE 


SHOCK WAVES 


Figure 1. Concial flow particle trajectories and stream surface (from Ref* 
1), and sketch of flow about an elliptic cone showing bow shock 
wave, cross-flow sonic surface and cross-flow shock waves. 










Figure 4* Computational boundaries and various possibilities for line 


relaxation. 
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Figure 5* The 64 x 64 mesh for a 18*39^t 3*17^ elliptic cone at 10^ angle of 


attack 
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^ PRESENT RESULTS 
O JONES (REF. 3) 



Figure 6b* Bow shock position for a curcular cone of 10^ half angle at 10^ 
angle of attack* Calculated using a 64 x 64 grid* Comparison is 

Q 

with the rotational calculations of Jones * 
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Flgure 7. Pressure variation between the bow shock wave and the body for a 
circular cone of 10® half angle at 10® angle with = 2.0. 

9 

Comparison is with the rotational calculations of Jones . 



Figure 8 


• Comparisons of the streamline patterns on a circular cone at 
a) 20^ and b) 10^ angle of attack with = 2.0. Note the lift- 
off of the leeward node as well as the formation of a supersonic 


zone in the cross flow. 
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A PRESENT RESULTS 
O SICLARI (REF. 4) 
GROSSMAN (REF. 5) 


AD o 


Figure 10. Comparison of the results using the Euler equations due to 
Siclarl^®, and a quasilinear formulation of the potential equation 
due to Grossman'^, with the present result for 
a 20^: 1.5® elliptical cone at 10® angle of attack with 
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